[1] 12 15 10 13 9 19 16 14 16 13 8 15 18 11 17 12 15 18 12 7 13 12 14 14 8
[26] 8 12 10 16 15 17 9 10 14 12 9
Rodney Dyer, PhD
A Raster is simply a matrix of values with some additional decorations on it that allow it to have a spatial context.
[1] 12 15 10 13 9 19 16 14 16 13 8 15 18 11 17 12 15 18 12 7 13 12 14 14 8
[26] 8 12 10 16 15 17 9 10 14 12 9
For each value in the matrix, when it is turned into a raster object:
The cell (pixel) has a defined spatial extent (width, height, & origin).
All the physical space represented by that cell has exactly the same value
The courseness of the raster is question dependent:
3x5 matrix for Continental US may not adequately capture elevation trends.
1m2 matrix for elevation may be a bit too big.
Notice that when I plot it out, it does not show the data, but a summary of the data along with some key data about the contents, including:
crs (missing)memory if the raster is not that big or out of memory if it is just referencing.By far, you will most commonly working with pre-existing raster data.
url <- "https://github.com/DyerlabTeaching/Raster-Data/raw/main/data/alt_22.tif"
r <- raster( url )
rclass : RasterLayer
dimensions : 3600, 3600, 12960000 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -120, -90, 0, 30 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : alt_22.tif
names : alt_22
values : -202, 5469 (min, max)
Notice that this raster has a defined CRS and as such it is projected and the extent relates to the units of the datum (e.g., from -120 to -90 degrees longitude and 0 to 30 degrees latitude).
Just like all things in R, raster objects can be visualized using built-in functions as well as functions from external libraries.
This particular raster is quite large (in terms of the number of cells)
class : RasterLayer
dimensions : 3600, 3600, 12960000 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -120, -90, 0, 30 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : alt_22.tif
names : alt_22
values : -202, 5469 (min, max)
These data only represent the elevation of the land. Where there is water, the value in the underlying matrix is NA.
| Cell Type | Count |
|---|---|
| Land | 2,469,350 |
| Water | 10,490,650 |
One of the first things to do is to crop the data down to represent the size and extent of our study area. If we over 10 million missing data points (the ocean) and most of Mexico in this raster above but we are only working with sites in Baja California (Norte y Sur), we would do well to excise (or crop) the raster to only include the area we are interested in working with.
extentextentLet’s use the beetle data from the Spatial Points Lecture as the data we will be working with.
library( sf )
library( tidyverse )
beetle_url <- "https://raw.githubusercontent.com/DyerlabTeaching/Raster-Data/main/data/AraptusDispersalBias.csv"
read_csv( beetle_url ) %>%
st_as_sf( coords=c("Longitude","Latitude"), crs=4326 ) -> beetles
beetles %>%
st_bbox() xmin ymin xmax ymax
-114.29353 23.28550 -109.32700 29.32541
Maybe rounding it to:
I’ll just use the old eyeball test to make the numbers ‘round’.
To crop the raster, we use the crop() function and it makes a new raster (and I can throw the old big one away).
class : RasterLayer
dimensions : 960, 840, 806400 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -116, -109, 22, 30 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : alt_22
values : -202, 2263 (min, max)
Notice the add=TRUE adds to the previous plot, and (2) Need to run whole chunk to see built-in plot overlays.
Masking is similar to cropping but with one main distinction. When you mask a raster, you do not reduce the size of it, you only allocate missing data to the part you are not interested in working with.
ggplotAs you probably guessed, there is a geom_raster() available to us. .redinline[However], we need to conver the data from a raster (matrix) to a data.frame object that ggplot can read.
raster to a data.frameA little coercion to move matrix into as.data.frame() is necessary. I also use the transmute() function which does in-place renaming (rather than mutate( X=y ) %>% select( -y ))
alt %>%
rasterToPoints() %>%
as.data.frame() %>%
transmute(Longitude=x,
Latitude=y,
Elevation=alt_22) -> alt.df
head( alt.df ) Longitude Latitude Elevation
1 -115.7958 29.99583 55
2 -115.7875 29.99583 126
3 -115.7792 29.99583 94
4 -115.7708 29.99583 99
5 -115.7625 29.99583 106
6 -115.7542 29.99583 120
geom_raster()You can work with raster data interactively (I just cannot do it here on this presentation because it has to be in real time).
click()Here are what the points look like.
Just like points, we can reproject the entire raster using the projectRaster function. Here I am going to project the raster into UTM Zone 12N, a common projection for this part of Mexico from epsg.io.
Unfortunately, the raster library does not use epsg codes so we’ll have to use the large description of that projection. See the page for this projection and scroll down to the proj.4 definition.
class : RasterLayer
dimensions : 981, 879, 862299 (nrow, ncol, ncell)
resolution : 834, 923 (x, y)
extent : -21583.09, 711502.9, 2428482, 3333945 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs
source : memory
names : alt_22
values : -202, 2222.994 (min, max)
What are the parts of Baja California that are within 100m of the elevation of site named San Francisquito (
sfran)?
To answer this, we have the following general outline of operations.
sfranalt raster that is within 100m (+/-) of that site.To do this we will use both the alt and the beetles data objects.
To extract data from raster objects, we need to coerce and specify.
library( ggrepel )
alt.df %>%
filter( Elevation >= 205,
Elevation <= 405 ) %>%
ggplot() +
geom_raster( aes( x = Longitude,
y = Latitude),
fill = "gray80",
data = alt.df ) +
geom_raster( aes( x = Longitude,
y = Latitude,
fill = Elevation) ) +
scale_fill_gradient2( low = "darkolivegreen",
mid = "yellow",
high = "brown",
midpoint = 305 ) +
geom_sf( aes(size=MFRatio),
alpha=0.5,
color="dodgerblue3",
data=beetles) +
geom_text_repel( aes( label = Site,
geometry = geometry),
data = beetles,
stat = "sf_coordinates",
size = 4,
color = "dodgerblue4") +
coord_sf()